Pacotes utilizados

library(plotly)
library(tidymodels)
## install BiocManager if not installed 
#if (!requireNamespace("BiocManager", quietly = TRUE))    
#install.packages("BiocManager")
## then install mixOmics
#BiocManager::install("mixOmics")
library(mixOmics)

Carregando a base de dados

Carregando base de dados para EHA/EHC

df <- read.csv('data2_areaR.csv',header = T, sep = ";",dec = ".")
rownames(df)= df[,1]                          
data <- df[,4:64]
categoria <- df$categoria
amostras <- df$amostra
bandeira <- df$bandeira

Análise exploratória dos dados

PCA

prin_comp <- prcomp(data, center = T, scale. = T)
components <- prin_comp[["x"]]
components <- data.frame(components, categoria, amostras)

#plotPCA
font.config <- list(
  family = "sans serif",
  size = 12,
  color = 'black')

fig <- plot_ly(components, 
        x = ~PC1, 
        y = ~PC2, 
        color = ~categoria, 
        colors = c('#BF382A', '#0C4B8E'), 
        type = 'scatter', 
        mode = 'markers',
        text = ~amostras)|> 
  layout(xaxis = list(title = "PC1 29,2%"), 
         yaxis = list(title = "PC2 17,8%"),
         title = "",
         font = font.config)
fig 
#orca(fig, file="PCA.png", scale = 3)
#orca(fig, file="PCA.svg")
tot_explained_variance_ratio <- summary(prin_comp)[["importance"]]
tot_explained_variance_ratio
##                             PC1      PC2      PC3      PC4      PC5      PC6
## Standard deviation     4.218463 3.294242 2.836302 1.905387 1.652718 1.416158
## Proportion of Variance 0.291730 0.177900 0.131880 0.059520 0.044780 0.032880
## Cumulative Proportion  0.291730 0.469630 0.601510 0.661030 0.705800 0.738680
##                            PC7      PC8      PC9     PC10     PC11      PC12
## Standard deviation     1.38593 1.243424 1.236495 1.095549 1.050944 0.9942217
## Proportion of Variance 0.03149 0.025350 0.025060 0.019680 0.018110 0.0162000
## Cumulative Proportion  0.77017 0.795520 0.820580 0.840260 0.858360 0.8745700
##                             PC13      PC14      PC15     PC16      PC17
## Standard deviation     0.9692859 0.9302106 0.8494444 0.814257 0.7787514
## Proportion of Variance 0.0154000 0.0141900 0.0118300 0.010870 0.0099400
## Cumulative Proportion  0.8899700 0.9041500 0.9159800 0.926850 0.9367900
##                             PC18      PC19      PC20      PC21      PC22
## Standard deviation     0.6944917 0.6296314 0.5886719 0.5500549 0.5179903
## Proportion of Variance 0.0079100 0.0065000 0.0056800 0.0049600 0.0044000
## Cumulative Proportion  0.9447000 0.9512000 0.9568800 0.9618400 0.9662400
##                             PC23      PC24      PC25      PC26      PC27
## Standard deviation     0.4718777 0.4624003 0.4404313 0.4373824 0.3995965
## Proportion of Variance 0.0036500 0.0035100 0.0031800 0.0031400 0.0026200
## Cumulative Proportion  0.9698900 0.9733900 0.9765700 0.9797100 0.9823300
##                             PC28      PC29      PC30      PC31      PC32
## Standard deviation     0.3620183 0.3287083 0.3105873 0.3031722 0.2879618
## Proportion of Variance 0.0021500 0.0017700 0.0015800 0.0015100 0.0013600
## Cumulative Proportion  0.9844800 0.9862500 0.9878300 0.9893400 0.9907000
##                             PC33      PC34      PC35      PC36      PC37
## Standard deviation     0.2678198 0.2398007 0.2296439 0.2085803 0.2007504
## Proportion of Variance 0.0011800 0.0009400 0.0008600 0.0007100 0.0006600
## Cumulative Proportion  0.9918700 0.9928100 0.9936800 0.9943900 0.9950500
##                             PC38      PC39      PC40      PC41      PC42
## Standard deviation     0.1930121 0.1781942 0.1746784 0.1704367 0.1633052
## Proportion of Variance 0.0006100 0.0005200 0.0005000 0.0004800 0.0004400
## Cumulative Proportion  0.9956600 0.9961800 0.9966800 0.9971600 0.9976000
##                             PC43      PC44      PC45     PC46      PC47
## Standard deviation     0.1526067 0.1441266 0.1319377 0.130614 0.1289309
## Proportion of Variance 0.0003800 0.0003400 0.0002900 0.000280 0.0002700
## Cumulative Proportion  0.9979800 0.9983200 0.9986000 0.998880 0.9991600
##                             PC48      PC49       PC50      PC51       PC52
## Standard deviation     0.1204568 0.1072623 0.08702224 0.0731475 0.07191255
## Proportion of Variance 0.0002400 0.0001900 0.00012000 0.0000900 0.00008000
## Cumulative Proportion  0.9993900 0.9995800 0.99971000 0.9998000 0.99988000
##                              PC53      PC54      PC55       PC56       PC57
## Standard deviation     0.06117772 0.0365507 0.0263285 0.02460594 0.02035427
## Proportion of Variance 0.00006000 0.0000200 0.0000100 0.00001000 0.00001000
## Cumulative Proportion  0.99994000 0.9999600 0.9999700 0.99998000 0.99999000
##                              PC58       PC59        PC60        PC61
## Standard deviation     0.01637325 0.01056493 0.009337684 0.008251142
## Proportion of Variance 0.00000000 0.00000000 0.000000000 0.000000000
## Cumulative Proportion  1.00000000 1.00000000 1.000000000 1.000000000
tit = 'PCA 3d'

fig1_3d <-plot_ly(components, 
        x = ~PC1,
        y = ~PC2,
        z = ~PC3,
        color = ~categoria,
        colors = c('#BF382A', '#0C4B8E'),
         text = ~amostras,
        type = "scatter3d",
        mode = "markers") |> 
  layout(scene = list(xaxis = list(title = "PC1 29,2%"), 
                        yaxis = list(title = "PC2 17,8%"),
                         zaxis = list(title= "PC3 13,2%")),
            title = "PCA")
fig1_3d
# orca(fig1_3d, file="PCA3d.svg")
fig2 <- plot_ly(components, 
        x = ~PC1, 
        y = ~PC2, 
        color = ~bandeira, 
        colors = c('orange','black','#BF382A','green4','#0C4B8E'), 
        type = 'scatter', 
        mode = 'markers',
        text = ~amostras)|> 
  layout(xaxis = list(title = "PC1 29,2%"), 
         yaxis = list(title = "PC2 17,8%"),
         title = "",
         font = font.config)
fig2 
# orca(fig2, file="PCA_band.png", scale = 3)
plot_ly(components, 
        x = ~PC1,
        y = ~PC2,
        z = ~PC3,
        color = ~bandeira,
        colors = c('orange','black','#BF382A','green4','#0C4B8E'),
         text = ~amostras,
        type = "scatter3d",
        mode = "markers")%>%
  layout(scene = list(xaxis = list(title = "PC1 29,2%"), 
                        yaxis = list(title = "PC2 17,8%"),
                         zaxis = list(title= "PC3 13,2%")),
            title = "PCA Bandeira")

PLS-DA

PLS-DA

Aqui, foi calculado uma PLS-DA levando em consideração a categoria para identificar qual são os marcadores importantes que distinguem cada marca.

X <- data
Y <- categoria

#Calculando a PLS-DA
result.plsda <- plsda(X, Y, scale = TRUE, ncomp = 5)

# CV
set.seed(123)
vald_plsda <- perf(result.plsda, validation = "Mfold", folds = 5, progressBar = F, auc = T, nrepeat = 100)
plot(vald_plsda, col = color.mixo(1:3), sd = T,legend.position="horizontal")

# https://rdrr.io/cran/mixOmics/man/plot.perf.html para ajudar no plot

# Escolhendo a melhor VL para a discriminação dos grupos
vald_plsda$choice.ncomp
##         max.dist centroids.dist mahalanobis.dist
## overall        4              4                4
## BER            4              4                4
# AUC
vald_plsda$auc$comp4
##    AUC.mean      AUC.sd 
## 0.968573000 0.003692385
# Plots
plsda_varcomps <- result.plsda$variates$X[,1:3]

datapls <- data.frame(plsda_varcomps,categoria)
cores <- c('#BF382A', '#0C4B8E')

# Plot interativo 2d
plot_ly(data = datapls,
        x = ~ comp1,
        y = ~ comp2,
        colors = cores,
        type = "scatter",
        mode = "markers",
        color = ~ categoria,
        size = 1)%>%
   layout(xaxis = list(title = "Componente 1"), 
         yaxis = list(title = "Componente 2"),
         title = "PLS-DA")
# 3d
plot_ly(data = datapls,
        x = ~ comp1,
        y = ~ comp2,
        z = ~ comp3,
        colors = cores,
        type = "scatter3d",
        mode = "markers",
        color = ~ categoria) %>%
  layout(scene = list(xaxis = list(title = "Componente 1"), 
                        yaxis = list(title = "Componente 2"),
                         zaxis = list(title= "Componente 3")),
            title = "PLS-DA 3d",
         font = font.config)
# Calculo VIPScores PLS-DA
vipscores <- vip(result.plsda) #calculo dos VIPScores da plsda 
vipscores.df <- data.frame(vipscores,colnames(X)) #df juntando vipscores + nome das variáveis
vip_cutoff = 1.1 #parâmetro

# Calculo dos VIPScores mais importantes para quarta componente
vipscores4comp <- vipscores.df[which(vipscores.df[,4] >= vip_cutoff),c(4,6)] #esse comando seleciona com o which quais variáveis são maiores que o parâmetro '=>vip_cutoff = 1'. o 'c(1,4) significa, o 1 a primeira componentes e o 4 a coluna com o nome da variável. Foi escolhido a 4 comp no comando 'vald_plsda$choice.ncomp'.

# Lista de Variáveis mais importantes
vipscores4comp
##          comp4 colnames.X.
## VOC10 1.247602       VOC10
## VOC11 1.307620       VOC11
## VOC13 1.190616       VOC13
## VOC25 2.837054       VOC25
## VOC28 1.170878       VOC28
## VOC29 1.175418       VOC29
## VOC31 1.158365       VOC31
## VOC32 1.166913       VOC32
## VOC33 1.171270       VOC33
## VOC34 2.054230       VOC34
## VOC36 1.110257       VOC36
## VOC37 1.190100       VOC37
## VOC38 1.182578       VOC38
## VOC40 1.216729       VOC40
## VOC41 1.210116       VOC41
## VOC42 1.195901       VOC42
## VOC44 1.695066       VOC44
## VOC59 2.184759       VOC59

A PLS-DA apresentou acurácia de 0,96857 com desvio padrão de 0,00369 na quarta variável latente. Excelente modelo para predição de EHA/EHC.

Foram selecionados na 4 componente, 27 variáveis com o VIPScore >= 1.

PLS-DA sem o grupo Shell

#Retirando o grupo shell da base de dados
df2 <- df[df$bandeira!='EHA-S',]
X2 <- df2[,4:64]
categoria2 <- df2$categoria
Y2 <- categoria2

#PLS-DA sem shell categoria
result.plsda2 <- plsda(X2, Y2, scale = TRUE, ncomp=5)

# CV
set.seed(1234)
vald_plsda2 <- perf(result.plsda2, validation = "Mfold", folds = 5, progressBar = F, auc = T, nrepeat = 100)
plot(vald_plsda2, col = color.mixo(1:3), sd = T,legend.position="horizontal")

#Selecionando a VL mais importante
vald_plsda2$choice.ncomp
##         max.dist centroids.dist mahalanobis.dist
## overall        3              3                4
## BER            3              3                4
#AUC da 3 comp
vald_plsda$auc$comp3
##    AUC.mean      AUC.sd 
## 0.963487000 0.003705924
#Plot 3d da PLS-DA Sem shell
plsda2_varcomps <- result.plsda2$variates$X[,1:3]
datapls2 <- data.frame(plsda2_varcomps,categoria2)
plot_ly(data = datapls2,
        x = ~ comp1,
        y = ~ comp2,
        z = ~ comp3,
        colors = cores,
        type = "scatter3d",
        mode = "markers",
        color = ~categoria2) %>%
  layout(scene = list(xaxis = list(title = "Componente 1"), 
                        yaxis = list(title = "Componente 2"),
                         zaxis = list(title= "Componente 3")),
            title = "PLS-DA 3d",
         font = font.config)
#Calculo VIPScores PLS-DA sem Shell
vipscores2 <- vip(result.plsda2) #calculo dos VIPScores da plsda 
vipscores.df2 <- data.frame(vipscores2,colnames(X2)) #df juntando vipscores + nome das variáveis
vip_cutoff2 = 1.1

#Calculo dos VIPScores mais importantes para primeira componente
vipscores3comp <- vipscores.df2[which(vipscores.df2[,3] >= vip_cutoff2),c(3,6)] 
#Esse comando seleciona com o which quais variáveis são maiores que o parâmetro '=>vip_cutoff = 1'. o 'c(1,4) significa, o 1 a primeira componentes e o 4 a coluna com o nome da variável. 

vipscores3comp
##          comp3 colnames.X2.
## VOC10 1.266392        VOC10
## VOC11 1.289764        VOC11
## VOC13 1.176945        VOC13
## VOC23 1.141260        VOC23
## VOC25 2.974478        VOC25
## VOC34 2.301792        VOC34
## VOC44 1.973896        VOC44
## VOC48 1.163187        VOC48
## VOC51 1.116619        VOC51
## VOC52 1.278763        VOC52
## VOC59 2.305147        VOC59

Foi escolhida a terceira componente para a PLS-DA EHA/EHC sem o grupo Shell e 11 compostos foram selecionados com VIPScores >1.1 para essa componente.

Determinação dos Marcadores por origem da aditivação

Marcadores selecionados - PLS-DA sem o grupo shell

# heatmap com melhor as variáveis da PLS-DA sem o grupo Shell
df.heatmap <- df2[,vipscores3comp$colnames.X2.]
df.heatmap_scaled <- scale(df.heatmap,center = TRUE, scale = TRUE) 
bandeira2 <- df2$bandeira
df.heatmap.cat <- data.frame(bandeira2,df.heatmap_scaled)

variables <- colnames(df.heatmap) 
values <- as.matrix(df.heatmap.cat[,2:12])

# cores
vals <- unique(scales::rescale(c(volcano)))
o <- order(vals, decreasing = FALSE)
cols <- scales::col_numeric("reds", domain = NULL)(vals)
colz <- setNames(data.frame(vals[o], cols[o]), NULL)

# heatmap bandeira
fig3 <- plot_ly(x=variables,y=bandeira2, z = values, type = "heatmap",colorscale = colz) |> layout(font=font.config)
fig3

Marcadores selecionados - PLS-DA

# heatmap com melhor as variáveis da PLS-DA sem o grupo Shell
df2.heatmap <- df[,vipscores4comp$colnames.X.]
df2.heatmap_scaled <- scale(df2.heatmap,center = TRUE, scale = TRUE) 
df2.heatmap.cat <- data.frame(bandeira,df2.heatmap_scaled)

variables2 <- colnames(df2.heatmap) 
values2 <- as.matrix(df2.heatmap.cat[,2:19])

# cores
vals <- unique(scales::rescale(c(volcano)))
o <- order(vals, decreasing = FALSE)
cols <- scales::col_numeric("reds", domain = NULL)(vals)
colz <- setNames(data.frame(vals[o], cols[o]), NULL)

# heatmap bandeira
fig5 <- plot_ly(x=variables2,
                y=bandeira, 
                z = values2, 
                type = "heatmap",
                colorscale = colz) |> layout(font=font.config)
fig5
#orca(fig4, file="heatmap.png", scale = 3)

Teste de performance dos marcadores selecionados

Teste com 3 marcadores x25,x29,x59

df.marcadores <- data.frame(categoria,df[,c(28,32,61)])
# Supondo que df.marcadores já foi criado conforme seu código
df.marcadores$categoria <- ifelse(df.marcadores$categoria == "EHA", 1, 
                                  ifelse(df.marcadores$categoria == "EHC", -1, df.marcadores$categoria))
# Transformando a coluna 'categoria' de volta para numérico
df.marcadores$categoria <- as.numeric(as.character(df.marcadores$categoria))

# Verificando a transformação
str(df.marcadores$categoria)
##  num [1:197] 1 1 1 1 1 1 1 1 1 1 ...
#Regressão
lm_model <- linear_reg() %>%
  set_engine('lm') %>% 
  set_mode('regression') %>%
  fit(categoria ~ VOC25 * VOC29 * VOC59, data = df.marcadores)
summary(lm_model)
##              Length Class      Mode
## lvl           0     -none-     NULL
## spec          7     linear_reg list
## fit          12     lm         list
## preproc       1     -none-     list
## elapsed       1     -none-     list
## censor_probs  0     -none-     list
#Predição
df_predict <- lm_model %>% predict(df)
rownames(df_predict) <- rownames(df) 
## Warning: Setting row names on a tibble is deprecated.
df_predict
## # A tibble: 197 × 1
##    .pred
##  * <dbl>
##  1 1.15 
##  2 0.477
##  3 0.762
##  4 0.808
##  5 0.970
##  6 0.883
##  7 1.48 
##  8 0.444
##  9 0.909
## 10 1.15 
## # ℹ 187 more rows
summary(lm_model)
##              Length Class      Mode
## lvl           0     -none-     NULL
## spec          7     linear_reg list
## fit          12     lm         list
## preproc       1     -none-     list
## elapsed       1     -none-     list
## censor_probs  0     -none-     list
#Plot
fig6 <- plot_ly(df, y = ~ lm_model$fit$fitted.values, x = ~categoria,
        colors = categoria,
        type = 'scatter', 
        alpha = 0.2, 
        mode = 'markers', 
        name = 'Tips',
        text = rownames(df_predict)) %>%
  layout(xaxis = list(title = "Categoria"), 
         yaxis = list(title = "Valor Predito"),
         title = "")

# Cálculo da acurácia
# Sendo, o positivo a aditivação:
# VP <- Qtd de EHA dado como EHA (qtd EHA > 0)
# VN <- Qtd de EHC dado como EHC (qtd EHC < 0)
# FP <- Qtd de EHC dado como EHA (qtd EHC > 0)
# FN <- Qtd de EHA dado como EHC (qtd EHA < 0)
df_predict1 <- data.frame(df_predict$.pred,rownames(df_predict), categoria)
VP <- df_predict1[which(df_predict1$df_predict..pred>0 & df_predict1$categoria=="EHA"),]
nrow(VP)
## [1] 91
VN <- df_predict1[which(df_predict1$df_predict..pred<0 & df_predict1$categoria=="EHC"),]
nrow(VN)
## [1] 98
FP <- df_predict1[which(df_predict1$df_predict..pred>0 & df_predict1$categoria=="EHC"),]
nrow(FP)
## [1] 0
FN <- df_predict1[which(df_predict1$df_predict..pred<0 & df_predict1$categoria=="EHA"),]
nrow(FN)
## [1] 8
AUC = 100*(nrow(VP)+nrow(VN))/(nrow(VP)+nrow(VN)+nrow(FP)+nrow(FN))
AUC
## [1] 95.93909
ESP=100*(nrow(VN)/((nrow(VN)+nrow(FP))))
ESP
## [1] 100
SENS=100*(nrow(VP)/((nrow(VP)+nrow(FN))))
SENS
## [1] 91.91919
#orca(fig6, file="fitted_values.png", scale = 3)

Das 8 amostras EHA que não foram classificadas corretamente, 4 eram de bandeira branca. Para testar vou adicionar um marcador para as “outras bandeiras”

Teste com 4 marcadores x25,x29,x59,x34

df.marcadores2 <- data.frame(categoria,df[,c(28,32,61,37)])
df.marcadores2$categoria <- ifelse(df.marcadores2$categoria == "EHA", 1, 
                                  ifelse(df.marcadores2$categoria == "EHC", -1, df.marcadores2$categoria))

# Transformando a coluna 'categoria' de volta para numérico
df.marcadores2$categoria <- as.numeric(as.character(df.marcadores2$categoria))

# Verificando a transformação
str(df.marcadores2$categoria)
##  num [1:197] 1 1 1 1 1 1 1 1 1 1 ...
#Regressão
lm_model2 <- linear_reg() %>%
  set_engine('lm') %>% 
  set_mode('regression') %>%
  fit(categoria ~ VOC25 * VOC29 * VOC59 * VOC34, data = df.marcadores2)
summary(lm_model2)
##              Length Class      Mode
## lvl           0     -none-     NULL
## spec          7     linear_reg list
## fit          12     lm         list
## preproc       1     -none-     list
## elapsed       1     -none-     list
## censor_probs  0     -none-     list
#Predição
df_predict2 <- lm_model2 %>% predict(df)
rownames(df_predict2) <- rownames(df) 
## Warning: Setting row names on a tibble is deprecated.
df_predict2
## # A tibble: 197 × 1
##    .pred
##  * <dbl>
##  1 1.02 
##  2 0.526
##  3 0.963
##  4 0.794
##  5 1.02 
##  6 0.878
##  7 1.35 
##  8 0.543
##  9 1.05 
## 10 1.11 
## # ℹ 187 more rows
summary(lm_model2)
##              Length Class      Mode
## lvl           0     -none-     NULL
## spec          7     linear_reg list
## fit          12     lm         list
## preproc       1     -none-     list
## elapsed       1     -none-     list
## censor_probs  0     -none-     list
#Plot
plot_ly(df, y = ~ lm_model2$fit$fitted.values, x = ~categoria,
        colors = categoria,
        type = 'scatter', 
        alpha = 0.2, 
        mode = 'markers', 
        name = 'Tips',
        text = rownames(df_predict2)) %>%
  layout(xaxis = list(title = "Categoria"), 
         yaxis = list(title = "Valor Predito"),
         title = "PLS-DA")
# Cálculo da acurácia
# Sendo, o positivo a aditivação:
# VP <- Qtd de EHA dado como EHA (qtd EHA > 0)
# VN <- Qtd de EHC dado como EHC (qtd EHC < 0)
# FP <- Qtd de EHC dado como EHA (qtd EHC > 0)
# FN <- Qtd de EHA dado como EHC (qtd EHA < 0)
df_predict2_1 <- data.frame(df_predict2$.pred,rownames(df_predict2), categoria)
VP2 <- df_predict2_1[which(df_predict2_1$df_predict2..pred>0 & df_predict2_1$categoria=="EHA"),]
nrow(VP)
## [1] 91
VN2 <- df_predict1[which(df_predict2_1$df_predict2..pred<0 & df_predict2_1$categoria=="EHC"),]
nrow(VN)
## [1] 98
FP2 <- df_predict1[which(df_predict2_1$df_predict2..pred>0 & df_predict2_1$categoria=="EHC"),]
nrow(FP)
## [1] 0
FN2 <- df_predict1[which(df_predict2_1$df_predict2..pred<0 & df_predict2_1$categoria=="EHA"),]
nrow(FN)
## [1] 8
AUC2= 100*(nrow(VP2)+nrow(VN2))/(nrow(VP2)+nrow(VN2)+nrow(FP2)+nrow(FN2))
AUC2
## [1] 95.93909

Acurácia ficou a mesma.

Teste com 4 marcadores x25,x29,x59,x11

df.marcadores3 <- data.frame(categoria,df[,c(28,32,61,14)])
df.marcadores3$categoria <- ifelse(df.marcadores3$categoria == "EHA", 1, 
                                  ifelse(df.marcadores3$categoria == "EHC", -1, df.marcadores3$categoria))

# Transformando a coluna 'categoria' de volta para numérico
df.marcadores3$categoria <- as.numeric(as.character(df.marcadores3$categoria))

# Verificando a transformação
str(df.marcadores3$categoria)
##  num [1:197] 1 1 1 1 1 1 1 1 1 1 ...
#Regressão
lm_model3 <- linear_reg() %>%
  set_engine('lm') %>% 
  set_mode('regression') %>%
  fit(categoria ~ VOC25 * VOC29 * VOC59 * VOC11, data = df.marcadores3)
summary(lm_model3)
##              Length Class      Mode
## lvl           0     -none-     NULL
## spec          7     linear_reg list
## fit          12     lm         list
## preproc       1     -none-     list
## elapsed       1     -none-     list
## censor_probs  0     -none-     list
#Predição
df_predict3 <- lm_model3 %>% predict(df)
rownames(df_predict3) <- rownames(df) 
## Warning: Setting row names on a tibble is deprecated.
df_predict3
## # A tibble: 197 × 1
##    .pred
##  * <dbl>
##  1 1.10 
##  2 0.563
##  3 0.957
##  4 0.804
##  5 1.01 
##  6 0.898
##  7 1.33 
##  8 0.537
##  9 0.881
## 10 1.13 
## # ℹ 187 more rows
summary(lm_model3)
##              Length Class      Mode
## lvl           0     -none-     NULL
## spec          7     linear_reg list
## fit          12     lm         list
## preproc       1     -none-     list
## elapsed       1     -none-     list
## censor_probs  0     -none-     list
#Plot
plot_ly(df, y = ~ lm_model3$fit$fitted.values, x = ~categoria,
        colors = categoria,
        type = 'scatter', 
        alpha = 0.2, 
        mode = 'markers', 
        name = 'Tips',
        text = rownames(df_predict3)) %>% 
  layout(xaxis = list(title = "Categoria"), 
         yaxis = list(title = "Valor Predito"),
         title = "PLS-DA") 
# Cálculo da acurácia
# Sendo, o positivo a aditivação:
# VP <- Qtd de EHA dado como EHA (qtd EHA > 0)
# VN <- Qtd de EHC dado como EHC (qtd EHC < 0)
# FP <- Qtd de EHC dado como EHA (qtd EHC > 0)
# FN <- Qtd de EHA dado como EHC (qtd EHA < 0)
df_predict3_1 <- data.frame(df_predict3$.pred,rownames(df_predict3), categoria)
VP3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred>0 & df_predict3_1$categoria=="EHA"),]
nrow(VP3)
## [1] 91
VN3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred<0 & df_predict3_1$categoria=="EHC"),]
nrow(VN3)
## [1] 98
FP3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred>0 & df_predict3_1$categoria=="EHC"),]
nrow(FP3)
## [1] 0
FN3 <- df_predict3_1[which(df_predict3_1$df_predict3..pred<0 & df_predict3_1$categoria=="EHA"),]
nrow(FN3)
## [1] 8
AUC3= 100*(nrow(VP3)+nrow(VN3))/(nrow(VP3)+nrow(VN3)+nrow(FP3)+nrow(FN3))
AUC3
## [1] 95.93909

Foi escolhido o Modelo 1 com três marcadores!